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04 ' We present our first numerical results of axisymmetric magnetohydrodynamic simula- 

^ tions for neutrino-cooled accretion tori around rotating black holes in general relativity. We 

^ , consider tori of mass ~ O.1-O.4M0 around a black hole of mass M = AMq and spin a = 0- 

. 0.9M; such systems are candidates for the central engines of gamma-ray bursts (GRBs) 

' formed after the collapse of massive rotating stellar cores and the merger of a black hole 

and a neutron star. In this paper, we consider the short-term evolution of a torus for a 
duration of ~ 60 ms, focusing on short-hard GRBs. Simulations were performed with a 
plausible microphysical equation of state that takes into account neutronization, the nuclear 
statistical equilibrium of a gas of free nucleons and a-particles, black body radiation, and 
a relativistic Fermi gas (neutrinos, electrons, and positrons). Neutrino-emission processes, 
I . such as capture onto free nucleons, pair annihilation, plasmon decay, and nucleon- 

O ' nucleon bremsstrahlung are taken into account as cooling processes. Magnetic braking and 

^ ' the magnetorotational instability in the accretion tori play a role in angular momentum 

C/3 , redistribution, which causes turbulent motion, resultant shock heating, and mass accretion 

. onto the black hole. The mass accretion rate is found to be M, ~ 1-10Mq/s, and the shock 

heating increases the temperature to ~ 10^^ K. This results in a maximum neutrino emission 
rate of = several x 10^^ ergs/s and a conversion efhciency LyjMtf? on the order of a few 
percent for tori with mass Mt « O.1-O.4M0 and for moderately high black hole spins. These 
results are similar to previous results in which the phenomenological a- viscosity prescription 



> 



' with the a-parameter of Ov = 0.01-0.1 is used. It is also found that the neutrino luminosity 



can be enhanced by the black hole spin, in particular for large spins, i.e., a <: 0.75M; if the 
accretion flow is optically thin with respect to neutrinos, the conversion efficiency may be 
<; 10% for a ^ 0.9M. Angular momentum transport, and the resulting shock heating caused 
0^ ' by magnetic stress induce time-varying neutrino luminosity, which is a favorable property 

for explaining the variability of the luminosity curve of GRBs. 

o 
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§1. Introduction 



There is growing evidence that gamma-ray bursts (GRBs) occur at cosmological 
distances. Assuming the isotropic emission of gamma-rays, the estimated absolute 
luminosity of many events is ^ 10^^ ergs/s. If the effect of anisotropy, such as 
the collimation of emission, is taken into account, the luminosity is estimated to be 
~ lO^*' ergs/s. Furthermore, the duration of the emission is very short, ~ 0.01- 
1000 s, indicating that the source is composed of a compact object of stellar-mass 
size.^^ Thus, the high luminosity can be explained by considering the conversion of 
energy from gravitational energy to thermal energy during accretion processes onto 
the compact object. Because a rotating black hole is the most efficient generator, it 
is now widely believed that the central engines of GRBs are composed of stellar-mass 
rotating black holes and massive, compact, and hot accretion disks (or tori). 

Over the past decade, many groups have studied the properties of dense, hot 
accretion disks (or tori) around a black hole [e.g., Refs. 2)-9)]. An often employed 
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method of study is to derive a stationary solution of accretion flow around a black 
hole, taking into account the detailed microphysics, but neglecting the effects of 
the thickness of the disks. ■^^"^^ Such studies have clarified qualitative and semi- 
quantitative features of the accretion disks, which are found to be dense, hot, and 
geometrically thick with maximum density ~ 10^^ g/cm^, maximum temperature 
~ 10^^ K, and geometrical thickness ~ 100 km for a black hole of mass ~ S-IOM©. 
Because of the high density and the high temperature, neutrino emission plays a 
crucial role in the dissipation of thermal energy in such accretion flows, resulting in 
a high neutrino luminosity of ^ 10^^ erg/s. One important property of this dense, 
hot accretion flow is that the density is so high that a large fraction of the neutrinos 
may be trapped by the flow and fail to escape before being swallowed by the black 
hole; i.e., the accretion flow is neutrino-dominated accretion flow (NDAF).^)"^) 

Although the stationary solutions reveal qualitative properties, these studies are 
not suitable for understanding the time variability of accretion flows. In the last three 
years, numerical simulations of viscous and neutrino-cooled accretion tori have also 
been carried out with detailed microphysics by two groups.®-*'^) Setiawan, Ruffert, 
and Janka have performed three-dimensional (3D) pseudo-Newtonian simulations 
with a tabulated equation of state (EOS) derived by Lattimer and Swesty^^^ and 
with a neutrino cooling employing a leakage scheme. As a transport mechanism of 
angular momentum in the flow, they incorporated a phenomenological viscosity using 
the so-called a-prescription.^^^ A general relativistic gravitational field is mimicked 
using a pseudo-Newtonian potential. They find that the efficiency of the conversion 
to neutrinos (defined by the ratio of the neutrino luminosity to the rest-mass-energy 
accretion rate) reaches typically to a few percent with maximum value of ~ 10%. 
Lee, Ramirez-Ruiz, and Page performed simulations similar to that of Setiawan et al., 
but with the assumption of axial symmetry, with slightly simplified microphysics, 
and for a longer time scale of ^ 0.5 s. They reported results similar to those of 
Setiawan et al. 

These simulations have provided a variety of nonstationary properties which 
cannot be found from the study of stationary accretion disks. They have found that 
the mass accretion rate and neutrino luminosity exhibit a mild time variability. Lee 
et al. found that the neutrino luminosity decreases on a time scale of 10-1000 ms, 
depending on the viscous parameter. They also reported that the inhomogeneity of 
the lepton number induces turbulent motion through a convective instability. These 
properties can be found only by numerical simulation. This implies that numerical 
simulation is a better approach for the study of GRB accretion flows. 

Despite their success, there are still many issues that have not yet been clarified 
by simulations. First, the simulations carried out to this time take into account gen- 
eral relativistic effects very crudely, using a pseudo-potential prescription. Although 
this method may partly capture the nature of general relativistic gravity around black 
holes, there is no guarantee that this potential can provide quantitatively accurate 
values concerning properties of black holes. Moreover, special relativistic effects are 
not incorporated in this method. The speed of matter around black holes is often 
a significantly large fraction of the speed of light, reaching a high Lorentz factor of 
~ 10, and hence it is reasonable to expect that the matter motion in the vicinity 



Magnetohydrodynamics of neutrino- cooled accretion tori around black hole 3 

of a black hole cannot be followed accurately with the pseudo-potential approach. 
Second, they add a so-called a-viscosity in their equations of motion to take into ac- 
count effects of angular momentum transport and resultant viscous heating. In the 
a-viscosity prescription, one has to a priori determine a value for the a-parameter 
(which is hereafter denoted by a^). This parameter essentially fixes the mass ac- 
cretion rate, the viscous heating rate, and the resulting neutrino luminosity. This 
implies that one never find these important quantities correctly, due to the unknown 
value of the a-parameter. Most people in the astrophysics currently believe that the 
magnitude of the viscosity in accretion flows is effectively determined by the turbu- 
lent motion of fluid induced by magnetic stress. If this is indeed the case, it implies 
that the actual effective viscosity varies in time. To capture the nature of a realistic 
accretion flow, a magnetohydrodynamics (MHD) simulation is needed, instead of 
adding an a-viscosity. 

With the motivation described above, we have performed a general relativistic 
magnetohydrodynamics (GRMHD) simulation using a code recently developed, 
which has already been applied to study the evolution of magnetized differentially 
rotating neutron stars-*^^)' -"^^^ and magnetized stellar core coUapse.^^^ In the present 
work, this code is used for a fixed black hole geometry but with detailed microphysics; 
a realistic EOS and neutrino cooling. The assumption of the fixed geometry is valid 
because we consider systems for which the mass of the torus is much smaller than the 
mass of the black hole. In such systems, the self-gravity of the torus does not play 
an important role. In the present context, magnetic stress natTirally induces time- 
varying angular momentum transport and shock heating. The rest-mass accretion 
rate and shock heating rate are computed in a first-principles manner. Furthermore, 
general relativistic effects are taken into account in a strict manner. Thus, it is 
possible to clarify the dependence of the neutrino luminosity in the GRB accretion 
torus on the spin parameter of black holes for the first time in dynamical simulations. 
In addition, we show that the magnetic effects induce a large time variability of 
the mass accretion rate and neutrino luminosity, which is not found in simulations 
employing the a-viscosity. 

This paper is organized as follows. In §2, we describe the basic equations, EOS, 
and neutrino processes that we take into account. In §3, the initial conditions and set- 
up for the simulation are described. In §4, numerical results are presented, focusing 
particularly on accretion rates of the rest-mass, energy, and angular momentum 
onto the black hole horizon and on the neutrino luminosity. Varying the masses of 
the torus and the black hole spin systematically, we clarify the dependence of the 
accretion rates and neutrino luminosity on these parameters. Section 5 is devoted to 
a summary and discussion. Throughout this paper, we adopt the geometrical units in 
which G = c = 1 where G and c are the gravitational constant and the speed of light. 
k and 7i are the Boltzmann and Planck constants. Latin and Greek indices denote 
spatial components and spacetime components, respectively. 77^^ and Sij{= 5^^) are 
the flat spacetime metric (in the cylindrical coordinates) and the Kronecker delta, 
respectively. 
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§2. Procedures for the numerical simulation 

We numerically solved GRMHD equations in the fixed geometry of a rotating 
black hole, assuming that the ideal MHD relation holds; the conductivity is assumed 
to be infinite. The basic equations and numerical methods are the same as those 
used in our first paper/^^ except for the EOS and microphysics that we employ. In 
this section, after we concisely review our basic equations for GRMHD, a detailed 
description of the EOS and microphysics is given. 

2.1. Definition of variables 

The fundamental fiuid variables are p, the rest-mass density, e, the specific 
internal energy, P, the pressure, and n^, the four velocity. For convenience, we also 
define the weighted rest-mass density p*, three- vector v^, specific enthalpy h, and a 
Lorentz factor w as 

P* = puW-9/v, (2-1) 

dt u^' ^ ^ 

h = l + e + -, (2-3) 
P 

w = au\ (2-4) 

where g is the determinant of the spacetime metric gf_iu , V is the determinant of ry^i/ , 
and a is the lapse function equal to (—5**)"^''^ (see ^2.2p . 

In addition, we need to define the electron fraction and the temperature T 
because these are used as arguments for tabulated EOSs (see ^2.41 for details). Here, 
Ye is defined by 

Ye = 2-5 
P 

where Ue is the number density of electrons (n_) minus that of positrons (n+), and 
m„ is the atomic mass unit. 

The only fundamental variable in the ideal MHD is b'^, a four- vector of magnetic 
field, that is perpendicular to , i.e., h^u^ = 0. ft'^ is regarded as the magnetic field 
observed in a frame comoving with the fluid. In the ideal MHD, the electric field E^^ 
in the comoving frame is zero and thus, the electric current is not necessary to 
determine the evolution of the field variables. Using b^, the electromagnetic tensor 
Ft"" is defined by^^) 

p^.u ^ ^^uap^^^^^ (2-6) 

where e^uap is the Levi-Civita tensor. 

Prom b^^, we define the three-magnetic field as 

= ay^{u%' - b^u'), (2-7) 

where 7 = 7/?? and 7 is the determinant of the three-metric 'jij. is regarded as 
the magnetic field observed in an inertial frame. We note that B^ = 0, Bi = "i-ijB^ , 
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and 



^ and bi = ^^(Bi + B^UiUi). (2- 



V-9/r] Wy/^ 

Using the hydrodynamic and electromagnetic variables, the energy- momentum 
tensor is written as 

Tf,u = {ph + b^)u^u^ + (P + lb^)g^.u - bf^b^, (2-9) 

where 

Then, the projection of the energy-momentum tensor is defined by 

PR = T^uu^n^ = {ph + b^)w^ - Ptot - {ab'f, (2-11) 
Ji = -T^yn^'Y^ = {ph + b^)wui - ab\, (2-12) 
Sij = Tf^^J^^'y^j = {ph + b^)uiUj + PtotHj - hbj, (2-13) 

where is the unit timelike vector normal to a spatial hypersurface, i.e., ^^yn" = 0. 
Also we define the total pressure Ptot = P+Pmagj where Pmag = b'^/2 is the magnetic 
pressure. 

Using Ji, and Sij, the energy- momentum tensor is rewritten in the form 

T^u = PRn^Uu + Ja\ny + Ja\n^ + Sij'y^'y^^. (2-14) 

This form of the energy-momentum tensor is useful for deriving the basic equations 
for GRMHD presented in ^2.31 For the following treatment, we define 

5*0 = V^PH and Si = \f^Ji. (2T5) 

We determined the evolution of these variables together with p^,, 1^, and B'' in the 
numerical simulation (see §2.30 . 

2.2. Gravitational field 

The GRMHD equations are solved in a fixed stationary gravitational field of a 
Kerr black hole. Following Refs. 16) and 17), we choose the Kerr-Schild coordinates 
of the Kerr black hole, because they have no coordinate singularity on the event 
horizon, and hence, regular solutions can be obtained even on the event horizon and 
slightly inside it. 

In cylindrical coordinates {Tu,z,ip), the line element of the Kerr-Schild solution 



is22) 



ds'^ = -dt^ + dw"^ + w'^dip^ + dz^ 

2Mr^ / rvjdw — aw'^dio zdz \ , , 

+ —A ^ 9 9 - + + dt\ , 2-16 
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where r is the radial coordinate in the Boyer-Lindquist coordinates, which is the 
positive solution of the equation 

^4 _ (^2 ^ ^2 _ ^2)^2 _ ^2^2 ^ 0_ (2-17) 

Here, M and a are the mass and the spin parameter of the Kerr black hole. A physical 
singularity of ring shape is located at w = a and z = 0, and the event horizon is at 
r = M + yfW — = "Th- We note that the determinant of the spacetime metric 
Qy^y has the simple form 

g = -vo'^. (2-18) 

For the Kerr-Schild solution, the lapse function a, the shift vector and the 
three-metric 7.^ are 



_ IMr^w _ 2Mr^zw _ 2Mar^ 

^ a{f + 2Mr^y' f + 2Mr^' ' a{f + 2Mr3)' ^ ' 

_ 2Mr^w'^ _ 2Mar^w^ 2Mr^wz 



n/ 2Ma^r^w^\ 2Mar^zw^ 2Mrz^ 

where / = + a^z^ and o" = + a^. The determinant of 7jj is 



7w = (^1 + j , Tv^^ = , = 1 + — ^ — , (2-21) 



g 2Mr^ 
— = TO ( 1 H 



7 = -^ = ^'(l + ^). (2-22) 



2.3. GRMHD equations 

Hydrodynamic equations to be solved ar^ 

Vf.ipu'') = 0, (2-23) 
7,^V^n = -Q;.7/, (2-24) 
n-V^n = -Qf.n''. (2-25) 

The first, second, and third equations are the continuity, Euler, and energy equations, 
respectively. Here, denotes a four-dimensional vector associated with neutrino 
cooling which is defined below, and is the covariant derivative with respect to g^i, . 
In the following, the equations are described in the cylindrical coordinates {w, z, ip) 
with the assumption of axial symmetry, and we define Sy = S^/ro, Uy = u^p/w, and 
By = B^vj. 



*' The equation V^T*^ = —Qi, is derived from the conservation equation of the total energy- 
momentum tensor T^^ , which is defined by the sum of energy momenta of all the matter fields. The 
derivation and assumption for the derivation is described in Appendix A. 
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With the quantities defined above, Eqs. (l2-23p ~ (l2-25p are written 



1 



dtp* H d^ip^v'^w) + dzip^v'^) = 0, 

w 



(2-26) 



dtSA + -a. 



Sav' + Ptot5/ 



1 



(n*)2 



B'{Ba + UAB'ui) 



-SodAa + SkdAp^ - ^SikdAf'' + —6^ 

Z VJ 



dtSy H kO', 



2(0 ^ 



,t\2' 



B^iBy+UyB'u 



5, 



dtSo + -d, 



Syv' 



1 



J-^B^ [By + UyB'Ui 



Qyi 



(2-27) 



(2-28) 



TU 



B'mB' 



B'uiB^ 



(2-29) 



where we use g = —w^. The subscript A denotes vo or and z, j, fc, • • • are tz7, z or 
99. Also, Kij is the extrinsic curvature, which is calculated in the stationary space 
from 



1 

2a 



(2-30) 



where is the covariant derivative with respect to jij. Equations ()2-26p - (j2-29p are 
solved in the method described in Ref. 12); we employ a high-resolution scheme 
consisting of the central scheme proposed by Kurganov and Tadmor^^^ with third- 
order cell-reconstruction. 

In addition to the hydrodynamic equations, we solve the evolution equation for 



dYe 
dt 



(2.31) 



where 7e is the capture rate of the electron whose definition is given in ^2.4[ Using 
the continuity equation, Eq. ()2-3ip can be rewritten as 



y^.{pYeU^^ 



-PU 7e 



(2-32) 



or, more explicitly. 



dt{p*Ye) H d^ip^Yev'^-oj) + dz{p*YeV^ 



w 



-P*le 



(2-33) 



We numerically solved Eq. ()2-33p in the same manner as Eq. (|2-26p . 
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When the ideal MHD relation holds, the Maxwell equations for the axisymmetric 
system are 



dtB^ 



-d,{B^v' -B^v^), 



dtB' 



1 



wiB^v' - B'v^ 



dtBy = d^iB^yy - Byy^) + d^{B^vy - Byy' 



(2-34) 
(2-35) 
(2-36) 
(2-37) 



Equation (j2-34p is the no-monopoles constraint, and Eqs. (|2-35p - (j2-37p are induction 
equations, which are solved to obtain the evolution of the magnetic field in the same 
manner as in Ref. 12), using the constraint transport scheme^^^ with a second-order 
accurate interpolation. 

The validity of the numerical code for the GRMHD equations was verified with 
several test simulations. The numerical results for standard test problems in rel- 
ativistic MHD, including special relativistic magnetized shocks, general relativistic 
magnetized Bondi flow in stationary spacetime, and a long-term evolution for a 
magnetized disk in full general relativity are presented in Ref. 12). 

After solving for the evolution of p*, Si, Sq, and together with B"^, we have 
to determine the primitive variables, such as p, e, Ui, and (or w = avf"). The first 
step in this procedure is to derive the following equation from the definition of Si: 



w 



D'^ihwy^iB'^ + 2hw). (2-38) 



Here, and are determined from the variables {p^,Si,B^) as 



B' 



B^ 



p*Vi 



and D'^ 



plVi ' 



(2-39) 



and to obtain Eq. (j2-38p . we use the relation SiB"^ = p^hB^Ui. Equation ()2-38p is 
regarded as a function of hw and 'w~'^ for given data set of s'^, B^, and D^. The 
definition of Sq can also be regarded as a function of hw and w^'^: 



So 

— = hw 
P* 



P* 



1 r 
2 



B'^w^'^ + D\hw 



,-2 



(2-40) 



Thus, Eqs. (j2-38p and (|2-40p constitute simultaneous equations for hw and w for 
given values of p*, 5j, 5o, and Y^, at each grid point. 

In our previous works^^^"^^^ in which the I^-law or hybrid EOSs are used, a single 
algebraic equation for hw can be constituted from Eqs. (j2-38p and (j2-40p . Then, 
a Newton-Raphson-type method can be used to obtain hw. In this work, we use 
tabulated EOSs for obtaining P, e, and h (these are functions of p, Y^,, and T; see 
^2.41 for details). In such a case, it is not easy to apply the same method, because of 
the complexity of the EOS. For this reason, we adopt a different iteration method. 

At each time step, and Ye are determined from their evolution equations. 
Then, p should be computed from p*/{y/jw), but w is the variable determined by 
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solving Eqs. (I2-38P and (l2-40p . Thus, we guess a solution for w as a first step and 
calculate a trial solution for p. For the trial value of w, we chose the value of the 
previous time step. For the resulting values of p with Y^, Eq. (j2-40p can be regarded 
as an algebraic equation for T, regarding h and P as functions of T. Then, searching 
for the solution from the table, we determine the value of T and, subsequently, P, 
e, and h. 

In the iteration, we have to determine a new trial value of w. For this purpose, 
we derive an algebraic equation for ur"^ from Eqs. (j2-38p and (j2-40p : 



If P and h are regarded as given values, this equation is a 3rd-order algebraic equation 



Cardano formula. 

We repeat these two procedures until a convergent solution is obtained. In most 
cases, the solution is obtained within ~ 10 iterations. However, in some cases, the 
solution is not obtained. The problem that often arises is that the solution for T is 
not found in the EOS table. This happens when e accidentally decreases to a small 
value for which the corresponding value of T is absent in the EOS tableEll In such a 
case, we set T to a minimum value that we choose arbitrarily. In the present work, 
we choose the minimum value to be 10^^^/^ K (see ^2.4p . 

2.4. Equation of state 

Dense, hot accretion tori of mass ~ O.l-IM© around a black hole of mass 3- 
4Mq are likely outcomes formed after the gravitational collapse of a massive stellar 
gQPg28)-30) after the merger of a low-mass black hole and neutron star,^^-*"^^^ as 
indicated by numerical simulations. According to the numerical results, the maxi- 
mum density and temperature of the tori are likely to be ~ 10^^ g/cm^ and T ~ 10^^ 
K. This temperature is high enough to photo-dissociate heavy nuclei, and hence, the 
main components of the baryon should be free protons, free neutrons, and a-particles. 
Because of this high temperature, electrons are relativistic and electron-positron pair 
production is possible in a low-density region. Furthermore, the degeneracy of elec- 
trons is high, due to the high density; the chemical potential of the electrons pe is 
comparable to or larger than kT. Photons are strongly coupled to the charged par- 
ticles, and hence, are trapped and advect with the matter flow. Neutrinos are also 
trapped in the high-density region with p ^ 10^^ g/cm^, whereas, for the low-density 
region, they escape freely from the tori and this contributes to the cooling process. 

*-* In the realistic EOSs for high-density matter with p > 10^ g/cm^, e is determined primarily 
by the radiation pressure for the high-temperature case, whereas it is determined by the electron- 
degenerate pressure for the low-temperature case. This implies that e is proportional to T^^^ at 
high temperature, whereas it approaches a constant for T <^ /ie/fc, where fie denotes the chemical 
potential of electrons. Therefore, if e drops below this limiting constant, no solution for T is found. 




(2-41) 



for w 



-2 



Thus, to obtain new trial value of w, we solve this equation using the 
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We determine the EOS, considering the physical conditions mentioned above. 
Om' approach is similar to that of Ref. 9). We assume that the torus is composed of 
free protons, free neutrons, a-particles, electrons, positrons, neutrinos, and radiation. 
Then, the pressure is written 

P = Pe + Pg+Pr + Pu, (2-42) 

where Pe, Pg, Pr, and P^ are the pressures of electrons and positrons, of baryons, of 
the radiation, and of neutrinos. 

In the following, we assume that the temperature is high enough (T ^ lO^'' 
K) that the electrons and positrons are relativistic (i.e., kT ^ rrieC^, where m-e is 
the electron mass). Then, the total pressure of the electrons and positrons can be 
determined analytically 



{kTY 



127r2(7i c) 



(2-43) 



where r]e = ^le/kT^ The number densities of electrons and positrons, n_ and n_|_, 
are related so as to give an electron fraction Yg of 



pYe {kTy 

— n_ — nj — 



(2-44) 



niu 3tt'^(?l c)^ 

We note that the assumption kT ^ mgC^ is not always valid, but for the dense region 
of the torus in which we are interested, it holds approximately.^^ Equations ()2-43p 
and (|2-44p hold for arbitrary degeneracy as long as the temperature is high enough. 
The gas and radiation pressure are written 

Pg = P^lll^^ (2.45) 

ruu 4 

Pr = (2-46) 

where X^nc is the mass fraction of free nucleons and the radiation density constant 
{'K'^k^ /Ibifi cf). If the neutrinos of all species are in thermal equilibrium with the 
matter, the neutrino pressure is given by 

P, = larT\ (2-47) 
o 

whereas the pressure is negligible if they are transparent. In this paper, we follow 
the treatment of Ref. 9) and set 

P, = larT\l-e-^^), (2-48) 
o 

where Ti, is the averaged optical depth of the neutrinos. 



*■* We reintroduce G and c in i) l2.4l and i]2.5l to clarify the dimensional units. 
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In the present context, the ratio of the radiation pressure to the gas pressure is 

^ = 3.0 X W-^Tf,p^l (2-49) 

where Tn = T/IO^^ K, pi2 = g/cm^, and we assume that X^^c = 1- For the 

case that the electron is degenerate with r]e ~ 10, the ratio of the electron pressure 
to the gas pressure is 

g^l.2xl0ir,W(^)'. (2.50) 

Thus, in the main body of the torus with Tn = 0(1) and pi2 = 0(1), the gas 
pressure is the dominant pressure source if the electrons are mildly degenerate with 
rje ^ 5, whereas for the case that the electrons are strongly degenerate, the degenerate 
pressure is dominant. 

The specific internal energy is written 

e = ee + Eg + Er + Eu, (2-51) 

where each term is related to the pressure by 

p " 2p p p 

The pressure and specific internal energy should be provided for given values of 
p, T, and Y^.^ For this purpose, the chemical potential of electrons, pe, and the free 
nucleon fraction, X^ua have to be written as functions of (p, T, Yg)- Equation (j2-44p 
is used to obtain the functional form of the chemical potential. Then, for simplicity, 
we use a fitting formula^^^ for the mass fraction of the free nucleon, 

Xnuc = Mm(22AT^l%-^/^ exp(-8.2/Tio), l) , (2-53) 

where Tio = r/10^° K and pio = p/W^^ g/cm^. 

We tabulated an EOS table for the ranges pmin < P < 10^^ g/cm'^ and Tmin < 
T < 5 X 10" K. The values 

Pmin aud Tjnin are the minimum values that we arbi- 
trarily chose to maintain the stability of the numerical computation. Because the 
temperature should not be much smaller than mgC^/Zc, because of our assumption 
that the electrons are relativistic, we set the value of Tj-nin to 10^+2/3 K. The value 
of 

Pmin is specified in ^3.4[ 
2.5. Neutrino emission 

Several processes contribute to the emission of neutrinos. The most important 
ones in the present context are the electron and positron captures onto free nucleons: 

p + e~ ^ n + i^e, (2-54) 
n + e+^p + Pe. (2-55) 



*' In Ref. 9), the authors made the additional assumption that Ye is a function of p and T. 
However, we do not make this assumption, as it does not always hold in the optically thin region 
for dynamical systems. 
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These processes change the electron fraction. Assuming that the electrons and 
positrons are relativistic, the electron capture rate in the optically thin region for 
neutrinos is approximately given by (see Appendix B) 

/ kT \^ 

7eO = Kc 2 [F4-{r]e)Xp - F4+{r]e)Xn], (2-56) 

\ TTlgC J 

where Xp and Xn are the mass fractions of free protons and free neutrons, = 
In 2/103 035 s-i, and 

F4-ix) = 45.59X + — x3 + -x^ + 24^6"^ - ^e'^^ + ^r^e'^^^ , (2-57) 



5 V 32 243 

F4+(x) = 24(e-- - + ^e'^^) . (2-58) 

In the optically thick region for neutrinos, /9-equilibrium is assumed to be satisfied 
and the electron fraction does not change. In this paper, we simply multiply by 
e""^" in order to take into account the opacity effect (diffusion effect) and define the 
electron capture rate by 

7e = 7eoe-"^ (2-59) 

We note that the actual optical depth depends on the neutrino species and neutrino 
energy. In this paper, we ignore such a dependence. 

As a result of the electron and positron captures, neutrinos are emitted and they 
carry away thermal energy of the rate 

6 



T77; Jxi'A. ' "\ 

Qcap = i^c — 2 [F5-iVe)Xp + F5+(r?e)X„], (2-60) 



where 



F5_(x) = 236.65 + + + 

ODD 

-120(e-^-^e-'^ + y^e-3-), (2-61) 

F5+(x) = 120(e-- - + ^e'^^) . (2-62) 

In addition to electron and positron captures, we consider the annihilation of 
electron-positron pairs into thermal neutrinos, nucleon-nucleon bremsstrahlung, and 
plasmon decay. For the neutrino emissivity due to electron-positron pair annihila- 
tion, wc use the fitting formula derived by Itoh et al.^^) For the other two, we write 
the emissivity as [e.g., Refs. 5), 36), and 37)] 

Qs = 1.5 X lO^^Tffpl^ ergs/cm3/s, (2-63) 
Qpia = 1.5 X 1032rf,7,,e-> (2 + 27^ + 7^) ergs/cm^/s, (2-64) 

where -fp = 5.565 x 10"^ ^(^i^Ts^fVS.^^) Then, we define the total emissivity as 

Q = (Qcap + Qpair + Qs + Qpla)e"^'' • (2-65) 
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Because Q is the emissivity (cooling rate) measured in the fluid rest-frame, we define 
as38)'39) 

Q;. = Qu^. (2-66) 

With this treatment, the effect of anisotropic emission associated with the matter 
motion is taken into account EH On the other hand, the effects of radiation transfer 
and the heating due to the emitted neutrinos are ignored. 

The main source of opacity for neutrinos is scattering by free nucleons and a- 
particles. For these processes, the order of the cross section is ~ 10~^^r^^ cm^.^''-' 
Because the number density of free nucleons is ~ 10^^ pn cm~^, the mean free path 
of the neutrinos is roughly 

~ 10^ T^lVn cm, (2-67) 

which is ~ 20GM/c^ for black holes of mass AMq. In this work, we consider accretion 
tori of width and thickness ~ 20-30GM/c^ and maximum temperature ^ 10^^ K. 
Hence, we simply set the optical depth to 

Tu = Cpu, (2-68) 

where C is a constant that we set to 1 or 1/3 [i.e., = p/{10^^ g/cm^) or p/(3 x 
10^^ g/cm^)]. Previous works investigating dense accretion tori (e.g., Ref. 9)) shows 
that the neutrinos are optically thick when the density satisfies the condition pn ^ 
1. Thus, the present treatment can capture the neutrino-trapping effect, at least 
qualitatively. 

2.6. Diagnostics 

We monitor the total baryon rest-mass M*, angular momentum J, internal en- 
ergy -Einti rotational kinetic energy T^ot, and electromagnetic energy E'eMi which are 
defined by 



= y p^^d-'x, (2-69) 

J = j P*hu^Vvd^x, (2-70) 

^int = j P*e^d^x, (2-71) 

Trot = j ^p^^hQu^u^^d^x, (2-72) 

Sem = j T^ua^d^x, (2-73) 

where is the electromagnetic part of the energy-momentum tensor, H is the 
angular velocity defined by u''^ /u^, and all the integrations are performed outside the 
event horizon. For the definitions of ii^int, Trot, and -EeMi we follow Ref. 14). 



*-* As a result of the anisotropic emission, angular momentum is dissipated, but the magnitude 
of this effect is much smaller than the loss associated with the infall into the black hole. 
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In stationary axisymmetric spacetime, the following relations are derived from 
the conservation law of the energy-momentum tensor in the absence of neutrino 
coohng: 



da 



(2-74) 
(2-75) 



From these relations, it is natural to define the energy and angular momentum 
accretion rates by the surface integral at the event horizon as 



E 
j - 



r=r-a 



r 



-gdS, 
/^dS, 



(2-76) 
(2-77) 



r=rii 



where dS = dOdtp. Prom the continuity equation, the rest-mass accretion rate is 
defined in the same manner as 



p^^v^r^dS. 



(2-78) 



Because we adopt cylindrical coordinates in the Kerr-Schild coordinates, grid points 
are not located at r = th- To obtain the values there, we use linear interpolations 
for all the necessary quantities. 

In the presence of neutrino cooling, Eq. (j2-74p is written 



~9n 



-gQut. 



Thus, the rate of energy loss by the neutrinos is given by 



-gutQd X. 



(2-79) 



(2-80) 



r>ru 



We refer to Lj^ simply as the neutrino luminosity. Strictly speaking, L^, is the sum 
of the energy emission rates of neutrinos toward infinity and toward the black hole 
horizon. Thus, the luminosity observed at infinity is smaller than this value, because 
a finite fraction of neutrinos emitted at each radius is always swallowed by the black 
hole. To roughly infer how much neutrino energy is likely to be swallowed by black 
holes, we also compute the quantity 



-gutQd'^x. 



where rph denotes the radius of the limiting circular photon-orbit 

•2 



2M 



1 + COS <! -cos"^(-a/M)| 



(2-81) 



(2-82) 



This radius is characterized by the fact that 50% of massless particles from a sta- 
tionary isotropic emitter at r = is captured by the black hole.^''^ To compute 
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the luminosity more accurately, it is necessary to multiply by an escape probability 
that depends on the position and velocity of the emitter. In the present paper, we 
do not take into account such probability but simply calculate the luminosity using 
Eq. (12^ or (|2^ . 

From M* and Lj,, the total rest-mass swallowed by a black hole and the emitted 
neutrino energy are defined by 

AM^{t) = [ ihdt\ (2-83) 

AE^,{t) = [ Lydt'. (2-84) 

In the following, we refer to AEy / AM^ as the average energy conversion efficiency 
to neutrinos. 



§3. Initial conditions and setting for computation 



3.1. Initial condition for torus 

As the initial conditions, we first prepare equilibrium states of a torus rotating 
around a black hole with no magnetic field and with no neutrino cooling. Such 
equilibria are determined from the first integral of the Euler equation, which has 
already been derived for axisymmetric stationary spacetime.'^'^^'^^^ We adopt the 
prescription of Ref. 41) because of its simplicity. In Kerr-Schild coordinates, the 
derivation of the basic equations is slightly more complicated than in Boyer-Lindquist 
coordinates. For example, in Boyer-Lindquist coordinates, we have = = 0, 
whereas they are nonzero in Kerr-Schild coordinates. We have to be careful regarding 
this point in setting the initial conditions. For completeness, here, we describe the 
basic equations for obtaining equilibria. 

First, we assume that the four- velocity has the following components: 

n'' = (n*,0,0,M^). (3-1) 

Note that this does not imply that u-^ and are zero, because of the presence of 
nonzero off-diagonal components of g^y. Defining the angular velocity i? = u'^/ti*, 
the Euler equation in the stationary axisymmetric system is written in the well- 
known form 

— 7 u'u^dkfi = 0. (3-2) 

w pa 

In the method of Ref. 41), one assumes 



£ = ^ = const. (3-3) 

Ut 



Then, Eq. (|3-2p can be rewritten as 



{dkUt)u\i-m)-^ = Q. (3-4) 

pri 
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Now, u^ut + u^u^ = —1 gives ti*(l — £f2) = —l/uf. Thus, we obtain 

dkH-ut) + ^ = 0. (3-5) 
ph 

If we assume that the fluid is isentropic, the first law of thermodynamics is used to 
rewrite the second term as dk ln{h). Thus, we obtain the first integral of the Euler 
equation in the same form as that in Ref. 41), 

hut = C, (3-6) 

where C is an integration constant. Then, using the relation 

^^_ 9ui + 9t^ ^ (3-7) 

we obtain 



Hence, h at each point is determined from Eq. (j3-6p for given values of £ and C. 

Specifically, the configuration of tori is determined by fixing the inner and outer 
edges on the equatorial plane. We denote the cyhndrical coordinate radii of these 
edges by wi and W2 (> roi). The condition h = 1 at w = wi and zu2 determines 
the values of C and £ using Eqs. (|3-6p and (|3-8p . Subsequently, the profile for h is 
determined from Eq. (|3-6p . 

To derive p and P from h, the EOS table is used. In doing so, we have to 
further assume relations among p, Yg, and T, because /i is a function of these three 
arguments in the table. In the present paper, we assume that T = 10^^ K uniformly 
in all regions. Because the tori and its atmosphere formed after stellar core collapse 
and the merger of compact objects are likely to be moderately hot at their formation, 
this assumption is reasonable. To determine Ye, we use the prescription of Ref. 9): 
For optically thick regions, we assume that /3-equilibrium holds, i.e., 

p + e~ ^ n + i/g, (3'9) 

and hence, the ratio of the proton fraction to the neutron fraction is assumed to be 

^ = eW-^^)A^, (3-10) 

where Q = 1.293 MeV is the difference between the neutron and proton masses. 
Here, we assume that the chemical potential of neutrinos is zero for simplicity. 

For optically thin regions, we assume that the reaction rates of the processes 
given in Eqs. (|2-54p and (j2-55p are identical and that the number densities of protons 
and neutrons are unchanged. '^^^ Under these assumptions, the electron fraction is 
approximately given by [see Eq. (12) of Ref. 9)] 



J^e — r Aj, 



1 Q/2 - 

- + 0.487- ' ^ 



2 kT 



(3-11) 
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where for X^^c, Eq. (j2-53p is used. In this rule with Tiq = 1, is smaller than 0.5 
for high-density regions with p ^ 10^^ g/cm^, whereas for low-density regions, it is 
larger than 0.5. 

The initial conditions obtained with the present method are not isentropic al- 
though we assume isentropic conditions in deriving Eq. (j3-6p . Thus, the torus is not 
exactly in equilibrium. However, in the test simulations with no magnetic field and 
no neutrino cooling, we found that the torus approximately remains in the initial 
state even after the density maximum orbits the black hole for 10 rotational periods. 
With this test, we confirmed that the initial conditions adopted in this method are 
approximately in equilibrium. 

There are several free parameters for determining the initial conditions: The 
black hole mass, M, the black hole spin, a, the mass of the torus, M^,, and the 
density profile of the torus. We fix the black hole mass to AMq. Such a black hole 
is a plausible outcome, forming soon after the stellar core collapse of massive star 
and after the merger of black hole-neutron star binaries. Because we do not know a 
plausible value of the black hole spin, numerical computations were systematically 
performed for a wide range of a. A plausible mass of a torus formed in the vicinity of 
a black hole is a few tenths of Mq. In the present work, we choose ~ O.1-O.4M0. 
Because we consider compact tori rotating around a black hole, the inner edge of the 
torus is chosen to be slightly outside of the radius of the innermost stable circular 
orbit (ISCO) as wi = 5-7M. Then, for a fixed value of M^,, the location of the 
outer edge is automatically determined to be ~ 30-40M for an I =constant angular 
momentum profile. 

In Table I, we list parameter values for selected models adopted in the nu- 
merical simulation. For models A-E, the torus mass is approximately given by 

~ O.25M0, whereas the black hole spins are all different. For models F~J, the 
black hole spin is fixed to a/M = 0.75, whereas the mass and inner edges of the 
torus are different. In Table II, we list the parameters for a test particle orbiting 
a black hole at the ISCO for several values of a/M. These are key quantities for 
determining the neutrino luminosity and the conversion efficiency of the rest-mass 
energy to neutrinos (see ^Ifoi' discussion). 

An important feature found from Table I is that for (approximately) fixed values 
of M, M*, and the radius of the torus, the maximum density is higher for larger values 
of a (compare models A-E). It is also found that for fixed values of M, M^,, and a, 
the maximum density is higher for smaller values of the rotation radius of the torus. 
For larger values of a, the radius of the ISCO on the equatorial plane is smaller 
(see Table II), and hence, the torus can be more compact, thereby increasing the 
maximum density. All these features imply that for larger values of a, the density can 
be higher. Because the optical depth of neutrinos depends strongly on the density 
for p ^ 10^^ g/cm^, this dependence of the maximum density on the black hole spin 
plays an important role in determining the neutrino luminosity. 

In Fig. [U^a), we plot density contour curves for model D (cf. Table I). It is seen 
that the torus has a geometrically thick structure in which the width ■072 — wi is 
comparable to the maximum thickness, ~ 30M. This is a universal property that 
holds for all the models. 
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Table I. Parameter set for the adopted models. The spin parameter of the black hole, the rest-mass 
of tori, the coordinate radius of the inner and outer edges, zui and zu2, of tori, the maximum 
density of tori, Eint/M,, Trot/M*, J/M^, and the rotation period at the density maximum, Pc, 
of tori. For every model, the black hole mass is fixed to be M = 4Mq, and the ratio of the 
electromagnetic energy to Ei^t is 1/200. zui and -072 are given in units of M, and -Eint and Trot 
are in units of M*/100. 



Model 


a/M 


M, (M(7,) 




/OinaxfE/cm"^) 
^iiitix Vo/ ^^^^ / 






J/M, 


Pc/M 




A 





0.244 


6.4, 41 


3.5 X 10'' 


0.85 


3.3 


4.03 


248 


1 


A2 





0.248 


6.4, 41 


3.6 X 10" 


0.85 


3.3 


4.03 


248 


1/3 


B 


0.25 


0.240 


6 , 36 


4.7 X 10" 


0.94 


3.7 


3.82 


219 


1 


B2 


0.25 


0.243 


6 , 36 


4.7 X 10" 


0.93 


3.7 


3.82 


219 


1/3 


C 


0.5 


0.249 


6 , 34 


5.4 X 10" 


0.99 


3.7 


3.70 


207 


1 


C2 


0.5 


0.253 


6 , 34 


5.5 X 10" 


0.99 


3.7 


3.69 


207 


1/3 


D 


0.75 


0.247 


6 , 32 


6.1 X 10" 


1.0 


3.8 


3.58 


203 


1 


D2 


0.75 


0.251 


6 , 32 


6.1 X 10" 


1.0 


3.8 


3.58 


203 


1/3 


E 


0.9 


0.247 


6 , 31 


6.4 X 10" 


1.1 


3.8 


3.51 


197 


1 


E2 


0.9 


0.251 


6 , 31 


6.5 X 10" 


1.1 


3.8 


3.51 


197 


1/3 


F 


0.75 


0.237 


4.8 , 24.6 


1.2 X 10'^ 


1.3 


4.9 


3.27 


146 


1 


G 


0.75 


0.132 


6.9 , 33 


2.8 X 10" 


0.82 


3.5 


3.74 


241 


1 


H 


0.75 


0.146 


6.6, 31.8 


3.4 X 10" 


0.87 


3.6 


3.68 


227 


1 


I 


0.75 


0.366 


5.4, 31.2 


1.0 X 10'^ 


1.2 


4.1 


3.45 


173 


1 


J 


0.75 


0.397 


6.0, 36 


7.8 X 10'^ 


1.1 


3.5 


3.61 


209 


1 



Table II. Quantities for a test particle orbiting a black hole at the ISCO for various values of 
the black hole spin. The value of the r coordinate of the ISCO, risco, the specific energy at 
the ISCO, -Eisco, and the specific angular momentum at ISCO, Jisco. The fifth and sixth 
columns list rn and rph. Note that the coordinate on the equatorial plane is related to r by 
zu = \/r^ + a'^. In this paper, we refer to 1 — Sisco as the specific gravitational binding energy 
at the ISCO. 



a/M 


nsco /Af 


Sisco 


Jisco /M 


rn/M 


rph/M 





6.000 


0.9428 


3.464 


2.000 


3.000 


0.25 


5.156 


0.9331 


3.221 


1.968 


2.695 


0.50 


4.233 


0.9179 


2.917 


1.866 


2.347 


0.75 


3.158 


0.8882 


2.489 


1.661 


1.916 


0.90 


2.321 


0.8442 


2.100 


1.436 


1.558 



Figure [2] plots the angular velocity Q as a function of the cylindrical radius on 
the equatorial plane. It is found that the torus has a differentially rotating velocity 
field with the steepness of the differential rotation \dln Q / din w\ ~ 2, steeper than 
the Keplerian law. This is also a feature of the i =constant velocity profile. 

3.2. Magnetic field 

To induce angular momentum transport during the evolution, a weak poloidal 
magnetic field is initially added to the equilibrium torus. In the present work, the 
profile of the magnetic field is chosen following Refs. 16)-18) as 



^10 tor p < pq. 
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X/M X/M 
(a) (b) 

Fig. 1. (a) Density contour curves of initial conditions for model D. The contour curves are drawn 
for p = ]^gi2-o.5i g^gjj^a with i = 1-10. The atmosphere is included in these initial conditions, 
(b) The same as (a) but for magnetic field lines. 




0.001 I ' — ' — ' ^ 

10 

X/M 

Fig. 2. Angular velocity as a function of the cylindrical radius on the equatorial plane for model 
D. 



where Aq is a constant that specifies the magnetic field strength and po is a cut-off 
density; the magnetic field is initially non-zero only in the region with p > pQ. Here, 
^0 is chosen such that the ratio of the magnetic energy to the thermal energy is 
1/200, and po is chosen as 0.2/>max) where pmax is the maximum density of the torus. 
Because the magnetic field is weak, the equilibrium configuration is not essentially 
modified by the magnetic force. 

In Fig. mb), we show the magnetic field lines for model D (cf. Table I). Com- 
paring Fig. m^a), it is seen that the magnetic field is initially present in a localized 
region near the density maximum because the density of the torus decreases rapidly 
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in the direction of the surface. 

Tori both with the poloidal magnetic field and with differential rotation are 
subject to winding of the field lines and subsequent magnetic braking. ^^-^ In addition, 
they are subject to the magnetorotational instability (MRI).^*^) Due to the winding 
of the magnetic field lines, the toroidal magnetic field grows linearly with time 
in the early phase of the numerical simulation as [e.g., Ref. 14)] 

B^ « tB^^ 2QtB^. (3-13) 

omw 

This amplification continues until the magnetic energy of the toroidal field becomes 
a few tenths of the rotational energy (cf. §4.1.ip . After the amplification stops, 
the magnetic braking transports the angular momentum outwards. For the present 
initial conditions with which -EEM/^rot ~ 0.001, the field growth saturates when B^ 
becomes ~ 2>QB'^ . This implies that the field growth stops at t ~ IbQ^^ ^ 2P^ 



O 5 



where Pq denotes the local rotational period. Because the magnetic field is confined 
near the density maximum, the expected saturation time is ~ 2-3 Po where Pc is the 
rotation period at the density maximum. 

The growth time (e-folding time) and wavelength of the fastest-growing mode 
of the MRI in the Newtonian theory are approximately given by [e.g., Refs. 24) and 
15)] 



^MRI = 2 



df2 



dlnw 



2'KV 

A 



4-1 -1/2 



(3-14) 
(3-15) 



where k is the epicyclic frequency of Newtonian theory. 



(3.16) 



and v\ = \jB^Bzl p'^ is the approximate Alfven velocity of the z component. In 
general relativity, the growth rate is modified, but it agrees with Eq. (|3-14p within 
~ 30%.^^) As in the case of the field winding, the growth time scale of the MRI is of 
order P^.. With the present initial conditions, we have w 10^2 ^nd P^ ~ 200M at 
the density maximum (see Table I), resulting in Amax ~ 2M. Because this is smaller 
than the width of tori, the MRI should play an important role in transporting the 
angular momentum and driving the turbulent motion in the accretion torus. We 
also note that the grid spacing is chosen to be typically 0.15M, so that the MRI 
can be resolved in numerical simulations. We note that with the present models as 
described above, we have -E'EMZ-E'int = 1/200. If the magnetic field strength is smaller 
than its present initial strength by a factor of more than 5-10, then the wavelength 
Amax is too Small to be resolved. If the magnetic field strength is chosen to be too 
large, Amax is so large that the MRI does not play an important role. The present 
choice of the field strength is appropriate for taking into account both the effects of 
the magnetic braking and the MRI. 
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It is known that non-magnetized tori with I =constant is often unstable with 
respect to the runaway instability (see, e.g., Ref. 43) and references cited therein). 
However, in the present model, the angular momentum of tori is redistributed by 
the magnetic braking and the MRI on a short time scale of ~ 2-3 Pc, and then a 
quasistationary accretion state, which has a less steep velocity profile, and hence, is 
stable with respect to the runaway instability, is established. The resulting quasista- 
tionary state is expected to depend only weakly on the initial profile of the angular 
momentum distribution, and thus, the numerical results should not depend strongly 
on the given velocity profile. 

3.3. Preparation of grid and boundary conditions 

A numerical simulation is performed in a nested grid zone; we divide the compu- 
tational domain into three zones and follow each zone with a different grid spacing. 
The computation in innermost zone is carried out with the finest grid spacing A. 
This zone covers the square domain satisfying < X < 48M, except for the inner 
square patch with < X < O.Srn, which is excised. Here, X denotes w and z. 
The second and third zones cover < X < 96M and < X < 192M, with the 
grid spacings 2A and 4Z\, respectively. At the outer boundaries of the first and 
second zones, values linear-interpolated from the 1-level larger zone are assigned as 
the boundary values. Specifically, the interpolation is performed for p^,, Ui = hui, 
T, Yg, and B^. Then, P and h are determined using the EOS table, and finally. Si 
and 5*0 are determined. At the outer boundaries of the third zone, we use an outflow 
boundary condition. The inner boundary condition at X = O.Srn, which is located 
well inside the event horizon, is arbitrary because it is inside the black hole, and 
backflow are automatically prohibited. 

The numerical simulation for the innermost zone is performed with A = 0.15M 
and 0.2M. For models D and E (see Table I), short-term simulations are performed 
with smaller grid spacings of Z\ = O.IM or 0.12M to check the convergence. We 
find that even for A = 0.2M, the torus is well resolved with the present model and 
parameter values. Thus, quantities such as the total rest-mass, the rotational kinetic 
energy, and the electromagnetic energy depend weakly on the grid resolution. The 
accretion rates at each time step, M^,, E, and J, which are evaluated at the event 
horizon of radius < 2M, depend more strongly on the grid resolution. However, 
an approximately convergent solution is obtained with A = 0.15M (see ^4.ip . The 
neutrino luminosity and the total emitted energy also depend strongly on the grid 
resolution, because the luminosity depends on a high power of the temperature, 
and hence a slight change in the temperature results in a significant change in the 
luminosity. However, we found that the total emitted energy converges within 20- 
30% for A = 0.15M. The convergence tends to be slower for larger values of a. The 
reason for this is that a larger value of a results in smaller radii of the event horizon 
and the ISCO, and thus a better grid resolution is required. 

3.4. Atmosphere 

Because any conservation scheme of hydrodynamics is unable to evolve a vacuum, 
we have to introduce an artificial atmosphere outside the tori. In particular, in the 
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MHD simulation, an artificial atmosphere of moderately large density is necessary 
because the computation often crashes when the ratio of the magnetic energy density 
to the rest-mass energy density exceeds a critical value. (In our present code, the 
critical value is ~ lO'^-lO^.) As the rule for introducing the atmosphere, we adopt 
the floor method proposed in Refs. 17) and 18). Specifically, when the density 
drops below a critical value, it is reset immediately. We choose a position-dependent 
critical density pcrit as 



Pcrit = Max 



p^ir.,Mm{p^i^{r/50M)-',300pmin) ■ (3-17) 



Our fiducial value for Pmin is 10^^/^ g/cm^. The floor-density prescription sacrifices 
the exact conservation of energy, rest-mass, and angular momentum. To elucidate 
the dependence of the numerical results on the value of the floor, we performed test 
simulations for model D using 

Pcrit = Max /Jmin,Min(^Pmm(?'/50M)"^, ISOpmin) , (3-18) 

with pmin = 10^°/^ g/cm^, and 

Pent = Max[/Jmin,Min(pmm(r/50M)-2^ lOOp^^in)] , (3-19) 

with Pmin = 10^^/^ g/cm^. We found that because the floor density is much smaller 
than that of the main body of the tori, the energy, rest-mass, and angular momentum 
are conserved to within ~ 10%. We also found that the total energy emitted by 
neutrinos and the total accreted rest-mass vary by no mor than ~ 10% as the floor 
values are varied. However, the density of atmosphere is still high enough to prevent 
formation of an outflow of a high Lorentz factor: During its outward propagation, 
the outflow carries a sufficiently large rest-mass and is thus decelerated. In other 
words, the velocity of the outflow depends strongly on the density of the atmosphere. 
In the present paper, we do not focus on the properties of the outflow. 



§4. Numerical results 



Numerical simulations were performed for a wide variety of models, listed in 
Table I, and for two or three grid resolutions in each case. The computations were 
continued to i ~ 60 ms (3000M). The total emitted neutrino energy and the total 
rest-mass swallowed by the black hole for < t < 50 ms are presented in Table III. 
As mentioned in ^3.41 these results depend on the density of the atmosphere weakly, 
varying by ~ 10%. 

Magnetized accretion torus evolves qualitatively similarly for a choice of the pa- 
rameter values. In ^4.1.H we summarize the general features and present the results 
for model D. In the subsequent subsections, we clarify the quantitative dependence 
of several quantities, such as the neutrino luminosity, maximum density, and maxi- 
mum temperature, on the mass and the initial rotation radius of the torus, the spin 
of the black hole, and the optical depth. 
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Table III. The total energy emitted by neutrinos in units of 10^^ ergs, the total accreted rest-mass, 
and the conversion efficiency. The time integration was performed for < t < 50 ms. The results 
with A/M — 0.2 and 0.15 are shown. For models with " — ", simulations were not performed. 
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4.1. General features: Results for model D 
4.1.1. Evolution of torus 

In Fig. [3l^a), we plot the evolution of the rotational kinetic energy and elec- 
tromagnetic energy for model D. For the first ~ 10 ms, magnetic field strength is 
amplified by the winding of the field lines and the MRI due to the presence of the 
poloidal magnetic fields and differential rotation. Here, the winding primarily con- 
tributes to the increase of the electromagnetic energy. This is found from the fact 
that the electromagnetic energy increases approximately in proportion to t^ for f ^ 5 
ms [see Eq. ()3-13p ]. The amplification continues until the electromagnetic energy 
reaches ~ 10% of the rotational kinetic energy. After the field growth saturates, 
transport of angular momentum and the resulting mass accretion onto the black 
hole are enhanced. The strong magnetic field subsequently induces turbulent matter 
motion driven by the MRI. As a result, shocks are generated and heat the matter, 
in particular for the inner part of the torus, up to ~ 10^^ K. This contributes to 
quick rise of the neutrino luminosity in the first ~ 10 ms (see Fig. [5]). At the same 
time, rotational kinetic energy is converted to thermal energy through the magnetic 
effects. Due to this effect, together with accretion, the rotational kinetic energy (as 
well as the internal energy) quickly decreases for t ^ 30 ms. For comparison, we 
show the evolution of the rotational kinetic energy obtained in the simulation of a 
non- magnetized torus. This reveals that the torus does not change much only in the 
presence of neutrino cooling. 

For t ^ 30 ms, no sharp decrease of the rotational kinetic energy nor sharp 
increase of the electromagnetic energy is observed. This indicates that the torus 
relaxes to a quasistationary accretion state. In this phase, the rotational kinetic 
energy gradually decreases as a result of the mass accretion. 
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Fig. 3. (a) Evolution of the rotational kinetic energy and electromagnetic energy in units of the 
initial value of the rest-mass, Mto, for model D with A/M = 0.1, 0.15, and 0.2. For comparison, 
the results with no magnetic field and with A/M = 0.2 are also shown, (b) The same as (a), 
but for the maximum density, (c) The same as (a) , but for the maximum temperature, (d) The 
same as (a), but for the maximum value of in units of the initial maximum value of B^. (e) 
The same as (d), but for the maximum value of . (f) The same as (a), but for the internal 
energy. For comparison, the electromagnetic energy for A — 0.1A4" is also shown. 



Figure El^b) plots the evolution of the maximum density. For 10 ms ^ t ^ 30 
ms, we see that it varies violently in the range between ~ 3 x 10^^ g/cm'^ and 
~ 2 X 10^^ g/cm^. This indicates that in such an early phase, the turbulent motion 
is strongly excited. For comparison, the result for the non-magnetized case is also 
shown. It is seen that the maximum density increases only gradually with time due 
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to neutrino cooling, showing that the torus is approximately in a quasistationary 
state. This indicates that the violent time variation of the maximum density in the 
early phase is induced by the magnetic effects. By contrast, the maximum density 
gradually increases for t ^ 30 ms. This is due to the fact that the rotation radius 
of the torus gradually decreases in the quasistationary mass accretion phase making 
the torus compact. [The density maximum is initially located at « WM, whereas 
it is at tz7 PS 5M at t = 40 ms; cf. Fig. EJa)]. 

Because the region with p ^ 10^^ g/cm^ is optically thick with respect to neutri- 
nos, a fraction of the neutrinos are trapped by the matter flow,^)'^) and this reduces 
the fraction of neutrinos that escape to infinity. This fact is clarified in ^4.41 which 
displays the dependence of the neutrino luminosity on The neutrino-trapping 
effect always plays a role for tori of mass ^ O.IMq in the model studied presently 
because pma.x of such tori is > 10^^ g/cm^ (cf. ^4.2p . 

Figure [St^c) plots the evolution of the maximum temperature. Because of shock 
heating, the temperature quickly rises in the first ~ 10 ms. Then, the maximum 
temperature reaches ~ 10^^ K and relaxes approximately to a constant value. This 
is a qualitatively universal feature for the tori considered in this paper. The value 
of the maximum temperature depends on the spin of the black hole and the mass of 
the torus, (cf. §M]and U^ . 

Figures [3jd) and (e) plot the evolution of the maximum values of \B^\ and {B^l 
in units of the initial maximum value of \B^\. In the early phase, the poloidal 
component of the magnetic field grows in an exponential manner, indicating that 
the MRI occurs. By contrast, the toroidal component {B^l increases linearly with 
time due to the winding of the magnetic field lines. This property is also seen in 
Fig. El^a), where the electromagnetic energy increases in proportion to t'^ for t ^ 5 
ms. After the amplification of the magnetic field, the electromagnetic energy density 
becomes comparable to the internal energy density [see Fig. [3)[|f)]. Then the growth 
of the field saturates. This is also one of the universal features for the evolution of 
magnetized tori, irrespective of the mass of the tori and the black hole spin. 

Figures [3l^b)-(e) show that numerical results for the maximum density, maximum 
temperature, and maximum magnetic field strength do not depend strongly on the 
grid resolutions. This indicates that the simulation with A = 0.15M, which is a 
typical choice for the grid resolution, provides a result of good convergence. We note 
that this is also the case for any value of a < 0.9M, although higher grid resolution 
is required for larger values of a, because the coordinate radius of the event horizon 
is smaller in this case. 

4.1.2. Evolution of accretion rates 

In Fig. m we plot the evolution of the rest-mass accretion rate, M^, accretion 
time scale, M^,{t) / M^^ and energy and angular momentum accretion rates, E and J, 
in units of M^. It is seen that as the electromagnetic energy increases [see Fig. [3^a)], 
the mass accretion rate quickly increases, and when the growth of the magnetic field 
saturates, it reaches ~ IOMq/s. Then, it gradually decreases and eventually relaxes 
to ~ I-2M0/S in a quasistationary state. We note that the accretion rate depends 
on the mass of torus and on the black hole spin (see ^4.21 and ^4.3p . 
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Fig. 4. Evolution of (a) the rest- mass accretion rate, M*, (b) the accretion time scale, A/*/M*, (c) 
the energy accretion rate in units of the rest- mass accretion rate, E/M*, and (d) the angular 
momentum accretion rate in units of the rest-mass accretion rate, J/M*, for model D with 
A/M = 0.1, 0.15, and 0.2. 



The associated accretion time scale for t ^ 30 ms is violently time-varying with 
the average value of ~ 50 ms, but it relaxes to ~ 0.1-0.2 s in the quasistationary 
state. In a previous work^^ in which the a-viscosity is used for the evolution of a 
neutrino-cooled torus of mass ~ O.SM©, the accretion time scale is ~ 50 ms for the 
a-viscous parameter value = 0.1 and ~ 0.5 s for = 0.01. The present numerical 
results approximately correspond to the case that 0.1 for the early phase and 

~ 0.04 for the quasistationary phase. This is a universal feature for any values of 
the mass of the torus and the black hole spin. 

The ratio of the energy accretion rate to the mass accretion rate, -E/M*, is 
~ 0.85-0.9 for the high-resolution runs. One may think that this is a reasonable 
result, because the ratio of the specific energy to the rest-mass energy of a test 
particle rotating around a Kerr black hole of a/M = 0.75 at the ISCO is 0.89 (see 
Table II). However, this is an accidental coincidence, as can be understood from the 
following two points, (i) The torus is not composed of test particles but fluid with 
a significant internal energy obtained by shock heating. Due to the contribution of 
the specific internal energy, the value of E should increase beyond that of the test 
particles, (ii) The accreting matter loses kinetic energy as it falls into the black hole 
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Fig. 5. Evolution of (a) the neutrino luminosity, L^, and (b) the efficiency of the conversion to 
neutrinos, Lu/M,, for model D with A/M = 0.1, 0.15, and 0.2. Here, "ph" denotes the result 
for the case that the luminosity is given by Eq. (l2-8ip . which approximately overlaps with Li,, 
as the difference is less than 10%. 



from the ISCO due to the magnetic stress. This effect should reduce the value of E. 
In the present case, these two effects approximately cancel each other. 

The ratio J /M^, is found to be ~ 1.5-2.3, which is smaller than the value of 
the specific angular momentum for a test particle orbiting at the ISCO, ~ 2.5 (see 
Table II). This indicates that during the infall of matter into the black hole from 
the ISCO, the magnetic stress extracts angular momentum, which is transported 
outwards along the field lines, as in the case of E. This effect was pointed out in 
Ref. 17). 

4.1.3. Neutrino luminosity 

Figure [5] plots neutrino luminosity L^, and the efficiency of the conversion to 
neutrinos (defined by the ratio of the neutrino luminosity to the rest-mass energy 
accretion rate, L^,/ Mj^c^). Soon after the magnetic field growth saturates at t ~ 10 
ms, the neutrino luminosity reaches a maximum of 2 x 10^^ ergs/s as a result of 
shock heating by violent matter motion induced by magnetic stress. Then, for 10 ms 
^ t ^ 30 ms, during which there is significant turbulent matter motion (cf. ^4.1.ip . 
it remains ^ 10^^ ergs/s. For t ^ 30 ms, the accretion rate gradually decreases, 
and so does the neutrino luminosity. In the quasistationary phase for t ^ 30 ms, we 
have Ly ~ 1-5 x 10^^ ergs/s. As shown in subsequent sections, these values for the 
neutrino luminosity depends on the mass of the torus and the black hole spin. 

We derive the neutrino luminosity by simply integrating the emissivity outside 
the event horizon. In an actual system of this kind, a fraction of the neutrinos 
emitted in the vicinity of black hole should be swallowed by it, because of the strong 
gravity. To estimate the dependence of the luminosity on the chosen domain for 
the integration, we also compute the luminosity, choosing the inner boundary of 
integration for the luminosity at r = rph (see the dotted curve in Fig. [5]). Note that 
only half of the neutrinos emitted by a stationary emitter at r = rph can escape. In 
this case, the luminosity is smaller by ~ 10% than L^,. Thus, the real luminosity 
may be ~ 90% of the value of Ly obtained from the integration for r > rn- 
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From the time integration, the total emitted neutrino energy, AE^, and the 
total accreted rest-mass for t < 50 ms are calculated, giving 3-4 x 10^^ ergs and 
~ O.1-O.13M0 for model D. (Note that these values depend on the grid resolution.) 
Thus, ~ 1.4-2.0% of the total accreted rest-mass energy, AMt,c^, is converted into 
neutrinos. (For A/M = 0.1, 0.15, and 0.2, it is 1.6%, 2.0%, and 1.4%, respectively). 
The order of magnitude of this value agrees with the results presented in Refs. 8) 
and 9). Our efficiency is smaller by a factor of ~ 2-3 than those obtained in the 
previous works. This is probably because of the differences in the treatments of 
the gravitational fields, neutrino opacity, and angular momentum transport process. 
Figure [5^b) also shows that the conversion efficiency is in a narrow range (between 
~ 1% and ~ 3%), close to the average value of Li^/M^c^ ~ 1-2%, for the entire time. 

For a/M = 0.75, the maximum hypothetical conversion efficiency is about 11% 
(cf. Table II). The results here imply that the conversion efficiency is not as large 
as the maximum. Note that the accretion time scale, M^/M^, ~ 100 ms, is approx- 
imately as long as the neutrino emission time scale, Ei^t/Ly. This implies that a 
part of the thermal energy is trapped by the matter flowing into the black hole and 
fails to be converted into neutrinos. Indeed, we find that the conversion efficiency 
depends strongly on the value of Q (see ^4.4p . This is one reason for the relatively 
small conversion efficiency. Another reason is that the initial radius of the torus 
is ~ lOM at the density maximum and at most 30M in the present model. The 
difference between the specific binding energy at r = lOM and at the ISCO is ~ 6%, 
and hence much smaller than 11%. We note that the radii of the ISCO are larger 
for smaller values of the black hole spin. Thus, the conversion efficiency is even 
suppressed for smaller values of the black hole spin, as shown in ^4.31 On the other 
hand, the suppression factor for larger spin with a/M = 0.9 is smaller. 

Numerical simulations indicate that the formation of a black hole-torus system 
from stellar core collapse and the merger of a black hole and neutron star is divided 
into two stages [e.g., Refs. 28)-33)]. In the first stage, the black hole and accretion 
torus are formed dynamically. According to our present numerical results, in such 
a stage, the amplification of magnetic fields and the subsequent redistribution of 
angular momentum proceed violently in the presence of magnetic fields of appreciable 
magnitude. As a result, shocks are generated, increasing the temperature of the 
torus, and the neutrino luminosity thereby reaches ^ 10^^ ergs/s. In the second 
stage, the system relaxes to a quasistationary state. In such a stage, the mass and 
spin of the black hole become approximately constant and the accretion proceeds 
with a time scale longer than that of the first stage. The present numerical results 
indicate that for this stage, the neutrino luminosity is of order 10^^ ergs/s. 

A characteristic feature found in the MHD simulation is that the luminosity 
curve is not smooth, in contrast to that reported in Refs. 8) and 9). This is due to 
the fact that in the present simulation, shock heating is associated with turbulent 
(i.e. irregular) matter motion driven by highly variable magnetic stress, whereas 
in the previous simulations, the heating is induced by a-viscosity, which leads to a 
nearly stationary accretion. 

The light curve of GRBs is often not smooth. If the GRBs are driven by pair 
annihilation of neutrino-antineutrino pairs emitted from the torus, the luminosity 
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curve of neutrinos should also be highly variable. The scenario in which the original 
thermal source of neutrinos is attributed to shock heating driven by the chaotic 
magnetic activity is thus favorable. 

The luminosity and total emitted neutrino energy, AE^, do not depend strongly 
on the grid resolution for a/M < 0.75. (The values of AEi, for A/M = 0.15 and 0.2 
agree within ~ 10-20%; cf. Table III.) This indicates that for computing the lumi- 
nosity, a grid resolution with A = 0.2M is acceptable for a < 0.75mE3 By contrast, 
for a = 0.9M, the luminosity and total energy for A/M = 0.15 and 0.2 do not agree 
well. (The relative error is by 20-30%; cf. Table III.) The reason for this difference is 
that the radius of the ISCO is much smaller than in the other models. (The ISCO is 
located at to ~ 2.5M; cf. Table II.) Because a large fraction of neutrinos are emitted 
near the ISCO, the luminosity depends strongly on the resolution near there. For 
a > 0.9M, a grid resolution with A = 0.2M is not acceptable for obtaining a result 
with good convergence. For a = 0.9M, we performed a simulation with A = 0.12M 
and found that the results with A = 0.15-M are in good agreement with those with 
A = 0.12M. This indicates that a grid resolution of A = 0.15M is acceptable even 
for a = 0.9M. 

4.1.4. Structure of the torus in the quasistationary phase 

In Fig. [6l we display color plots and contour curves for the density, temperature, 
electron fraction, neutrino emissivity, and ratio of the magnetic pressure to the gas 
pressure at t ~ 40 ms for model D. At t = 40 ms, the accretion torus relaxes to 
an approximately quasistationary state. The density maximum of the torus with 
Pmax ~ 10^^ g/cm^ is located at ~ 5M on the equatorial plane. Only a small 
part of the inner region with zu ^ lOM and with \z\ ^ 2M has a density larger 
than 10^^ g/cm^ and is optically thick. The temperature is also highest, ~ 10^^ K, 
for this high-density region, because the cooling is not efficient for such a region, 
due to its high opacity . The temperature of the optically thin region with p < 
10^^ g/cm^ is ~ 3 X 10^*^ K near the equatorial plane. An interesting feature is that 
the temperature is not uniformly low in the region above the torus. The reason 
for this is that the magnetic stress in the torus induces an outflow which ejects 
gas with high temperature. This dilute gas subsequently emits neutrinos, and the 
temperature thereby quickly decreases. 

As reported in Ref. 9), the electron fraction inside the torus is low, ^ 0.1. 
This reflects the fact that the density is so high that the electrons are highly de- 
generate. A typical value of the degeneracy parameter is rye ~ 2-4 in the region 
with p ^ 10^^ g/cm^. Around the envelope of the torus in the low density region 
{p ^ 10^ g/cm^), by contrast, the degeneracy is low, with rje < 1 for a large fraction 
of the fluid. 

Neutrinos are efficiently emitted in the region satisfying w ^ lOAf and \z\ ^ 5M. 
The neutrino emissivity is highest near the surface of the torus (not at the density 



*' The mass accretion and neutrino emission rates are induced by magnetic stress, which causes 
turbulent motion. As a result, the numerical results for different grid resolutions do not agree at 
each moment of time. However, the average values over a duration ~ 10 ms do not depend strongly 
on the grid resolution. 
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Fig. 6. Color plots for (a) the density, (b) the temperature, (c) the electron fraction, (d) the neutrino 
emissivity, and (e) the ratio of the magnetic pressure to the gas pressure at t ~ 40 ms for model 
D with A/M = 0.15. For (d), the region with Q < lO^'' ergs/s/cm^ is shaded black. The 
contour curves are plotted for (a) p = 10*"''' g/cm^ {i = 0-3), (b) log(rio) = 0, 0.4, and 0.8, (c) 
Ye = 0.1 X j (i = 1-5, (d) Q = 10^*+' ergs/s/cm^ (i = 0, 3), and (e) P^ag/P = iQ-^+^i (i = 0-2). 
(f) The angular velocity as a function of the cylindrical radius on the equatorial plane (solid 
curve). The dotted line denotes MQ = (M /■cu)'"''^'^ . Note that M ~ 5.9 km in the geometrical 
units. 
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maximum), because the region in the vicinity of the density maximum is so dense 
that the neutrinos cannot escape. In particular, the region near the inner edge of 
the torus has high emissivity. A large fraction of the neutrinos emitted from there 
will propagate toward the symmetry axis, where the neutrinos can escape freely, and 
hence, the pair annihilation rate of neutrinos and antineutorinos is expected to be 
highest in the vicinity and around the symmetry axis of a rotating black hole. 

The distribution of Pmag/P exhibits a clear contrast. In the torus, the value of 
Pma.g/P is much smaller than unity, implying that the gas pressure is dominant. By 
contrast, the value near the symmetry axis is much larger than unity, and hence, the 
gas pressure is negligible there. This structure has already been found in previous 
GRMHD simulations with a simple J^-law EOS.^'')'^^-' The present results show that 
the formation of such a distribution of Pm&g/P does not depend on the EOS and the 
neutrino cooling. 

Figure [6|^f) plots the angular velocity along the cylindrical radius in the equato- 
rial plane. This illustrates that the velocity of the torus is approximately Keplerian 
in the quasistationary phase. At the beginning of the simulation, the profile of the 
angular velocity is steeper than this, as shown in Fig. [2j Due to the magnetic 
braking and the MRI, the angular momentum is redistributed, and the profile of the 
angular velocity is modified and becomes Keplerian. The resulting torus is likely to 
be stable with respect to nonaxisymmetric instabilities and the runaway instability 
because of the Keplerian velocity profile. 

4.2. Dependence on the mass and the initial radius of the torus 

In Fig. [71[a), we plot the neutrino luminosity as a function of time for models D, 
G, and J. The mass of torus varies among these models, while the spin parameter of 
the black hole is the same and the initial rotation radii of the torus are approximately 
equal at t = 0. Figure [TU^a) shows that the luminosity systematically increases 
as the mass of the torus increases. This is simply due to the fact that there are 
more emitters for the larger-mass models. Because the maximum density is slightly 
higher (see Table I), there should be more trapped neutrinos for the larger-mass 
models, suppressing the neutrino luminosity. However, this does not decrease the 
neutrino luminosity significantly. The temperature is slightly higher for the larger- 
mass models [see Fig. El^b)]. This is one reason that the neutrino luminosity increases 
with the torus mass. 

The mass accretion rate is larger for the larger-mass models in the early phase, 
i.e., for t ^ 40 ms [see Fig. [7)^c)]. As a result, at late times {t ^ 40 ms), the differences 
among the mass accretion rates of the three models are relatively small, resulting 
in relatively small differences among the neutrino luminosities. The authors of Ref. 
8) report that the neutrino luminosity is approximately proportional to for the 
case of nonzero a-viscosity. In the present result, the maximum luminosity and total 
energy carried by neutrinos are approximately proportional to M^,. 

Figure E^d) plots the conversion efficiency, Ly/M^,, as a function of time for 
models D, G, and J. Although the neutrino luminosity differs significantly among 
three models, the conversion efficiency, varying between ~ 1% and ~ 5%, does not 
depend on the mass of the torus as strongly as the luminosity. The average conversion 
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Fig. 7. Evolution of (a) the neutrino luminosity, (b) the maximum temperature, (c) the mass 
accretion rate, and (d) the efficiency of the conversion to neutrinos L^/M*, for models D, G, 
and J with A/M = 0.15. Because the initial rotation radii are slightly larger, the luminosity 
and temperature for model G rise at slightly later time than for other models. 



efficiency, AEu/AM^,, is also in a narrow range between 1.4% and 1.8% for models 
G-J (see Table III). 
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Fig. 8. Evolution of (a) the mass accretion rate, M«, and (b) the neutrino luminosity, L^, for 
models D and F with A/M = 0.15. 



To see the dependence of the rest-mass accretion rate and neutrino luminosity 
on the initial radius of the torus, we performed a simulation for model F. This model 
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has the same black hole spin and approximately the same torus rest-mass as model 
D, but it has a smaller initial rotation radius. In Fig. [HI we plot the mass accretion 
rate and the neutrino luminosity for models D and F. Because model F has a smaller 
radius initially, the mass accretion rate and the neutrino luminosity quickly increase 
at earlier times. Also, the peak value of the luminosity for model F is larger by a 
factor of ~ 2 than that for model D, whereas the luminosity at late times is smaller 
for model F. However, the total emitted neutrino energy and the efficiency of the 
conversion to neutrinos are similar (see Table III). 

4.3. Dependence on the black hole spin 

Figure [9] plots the mass accretion rate and the neutrino luminosity as functions 
of time for models B-E. The black hole spin is different for each of these models, 
but the mass and rotation radii of the torus are approximately equal at t = 0. We 
find the following, (i) The luminosity for model E is largest among the four models, 
whereas the mass accretion rate is smallest for model E, probably because its event 
horizon has the smallest area, (ii) The luminosity increases with a by a significant 
amount only for a > 0.5M (see also Table III). The magnitude of the luminosity 
varies only slightly for a/M = 0-0.5, while the luminosities for models B and C are 
approximately the same, (iii) The decay time of the luminosity increases with a for 
a/M > 0.5. The facts (i)-(iii) indicate that for a sufficiently large value of the spin, 
a ^ 0.75M, the effect of the black hole spin enhances the neutrino luminosity, but 
for smaller values, the luminosity is not significantly enhanced by the spin effect. In 
the following, we describe the reasons for these types of behaviors. 

There are two factors which determine the dependence of the luminosity on 
the black hole spin. One is the fact that with larger spin, the radius of the ISCO 
decreases as the spin increases (cf. Table II). Then, the gravitational binding energy 
at the ISCO (1 — -^1300)5 which could be converted to thermal energy, increases as 
a function of a (cf. Table II). Indeed, the maximum temperature also increases as 
a increases [see Fig. [TOTa)]. This results in an enhancement of the efficiency for 
converting gravitational binding energy into neutrinos. Note that the gravitational 
binding energy at the ISCO changes slowly as a function of a for small values of a. 
This implies that this effect is not very important for a <C M. 

The small radii of the ISCO for larger values of the spin also result in the 
decrease of the mass accretion rate and in the increase of the accretion time. To 
achieve a high neutrino luminosity, a longer accretion time is favorable, because a 
fraction of the thermal energy trapped with the matter flowing into black hole is 
reduced. Thus, the case of a large spin has an advantage with regard to enhancing 
the neutrino luminosity due to (i) the large binding energy at the ISCO and (ii) the 
long accretion time. 

The other factor is that for a larger value of a, the rotation radius of the torus 
can be smaller. This results in a more compact torus which has a higher temperature 
and larger density. As a result of the larger density, (i) the trapped fraction of neu- 
trinos is increased, whereas (ii) the neutrino luminosity could be increased because 
of the higher temperature and larger density. Figure [10] plots the evolution of the 
maximum temperature and density for models B-E. For model B, the maximum 
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density gradually decreases below 10^^ g/cm^, which implies that no trapped neu- 
trinos exist in the late phase. For model C, pmax is slightly larger than 10^^ g/cm^, 
and hence, a small fraction of neutrinos are trapped. For models D and E, in which 
Pmax ~ 10^^ /cm^, a large fraction of neutrinos are trapped, whereas the emissivity 
is enhanced by the large density!*^ 

For models B and C, the above two factors seem to play an accidently identical 



role, resulting in similar luminosities !**^! For models D and E, on the other hand, 
the enhancement of the luminosity at high temperature and high density plays a 
stronger role than the effect of the opacity. Because the luminosity for models D 
and E is much larger than for models A-C, the longer accretion time and high 
temperature resulting from the large spin play the most significant role in enhancing 
the luminosity. 



*' The increase of the maximum density as a function of the spin is partly due to the fact that 
the mass accretion rate decreases as a function of the spin. 

**' In the discussion of this section, we use numerical results obtained with C = 1- K we instead 
used ^ = 1/3, the situation would be different. In this case, the fraction of trapped neutrinos is 
small for model C2, and hence, the luminosity for this model is larger than that for model B2. 
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In the present numerical work, there is the additional reason that the conversion 
efficiency increases rapidly with the spin for a/M ^ 0.75. As mentioned above, the 
radius of the ISCO for black holes of larger spins is smaller, and a larger fraction 
of the rest-mass energy can be converted into thermal energy (cf. Table II). In the 
present model, the density maximum for the torus is initially located at ~ lOM. 
For smaller values of the spin, then, the possible conversion rate for matter near the 
density maximum is much smaller than the maximum allowed value given in Table 
II, because the difference between the specific binding energy at the ISCO and at 
the density maximum is much smaller than the specific binding energy at the ISCO. 
For large values of the spin with a — > M, by contrast, the difference is close to the 
specific binding energy at the ISCO. This is one reason that the conversion efficiency 
for models D and E is much larger than those for other models. 

4.4. Dependence on ( 

In this subsection, we compare the results for pair models for which the black 
hole spin and properties of the torus are identical but the rule for identifying the 
optically thick region for neutrinos is different. We determine that the region with 
Pii > 1 is optically thick for models A-E, whereas the region with pn > 3 is for 
models A2-E2. 

Figure [TlTa) plots the evolution of the neutrino luminosity, Ljy, and the con- 
version efficiency, L^/M^, for models D and D2. Because the optically thin region 
is wider, the neutrino luminosity for model D2 is higher than that for model D, in 
particular in the early phase, with t ^ 30 ms. The relative difference is ~ 30-40%. 

The conversion efficiency, Lj^/M^,, for model D2 is larger than that for model D, 
reflecting the higher luminosity for model D2. It often reaches 10% and is ~ 2.5%(~ 
AE^/AM^:C^) on average. Nevertheless, it is smaller than the maximum value, ~ 
11%. Recalling that the accretion time scale, M^/M^ ~ 100 ms, is approximately as 
long as the neutrino emission time scale, Ei^t/L^, this implies that a large fraction 
of the neutrinos are still trapped by the matter and advected into the black hole, 
thus failing to escape from the torus. Indeed, the maximum density for model D2 is 
larger than 3 x 10^^ g/cm^ for most of the time [see Fig. [TTTd)]. 

Figure [TlTc) plots the evolution of the neutrino luminosity for models A and A2. 
As in models D and D2, the luminosity for model A2 is larger than that for model A. 
However, the difference is not as large as that between models D and D2. The reason 
for this smaller difference is that the density of the torus for models A and A2 is not 
as large as that for models D and D2 [see Figs. [TOlfb) and IllTd)]. In particular, for 
t ^ 20 ms, the density is smaller than 10^^ g/cm^, and hence, no region is optically 
thick for models A and A2. We conclude that for a torus with small values of a and 
with mass ~ O.25M0, the neutrino luminosity depends only weakly on the value of 
C 1/3). 

The same conclusion is reached for tori of small mass which would have densities 
smaller than 10^^ g/cm^. (Compare the maximum densities for models D, F-J in 
Table I.) Actually, only a small fraction of a torus with M^, ^ O.IMq is optically 
thick for a < 0.75M. For a ~ 0, the neutrino-trapping effect plays an important role 
only for massive tori with M* ^ 0.3Mq. For large values of a, i.e., for a ^ 0.75M, 
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Fig. 11. Evolution of (a) the neutrino luminosity, L^, and (b) the efficiency of the conversion to 
neutrinos, L^/M,, for models D and D2 with A/M — 0.15. (c) The same as (a), but for models 
A and A2 with A/AI = 0.2. (d) Evolution of the maximum density for models A2, D2, and E2. 
(e) The same as (a), but for models E and E2 with A/M — 0.15. (f) The same as (b) but for 
models E and E2 with A/M = 0.15. 



by contrast, neutrino trapping plays a role even for M^, ~ O.IMq. 

In Fig. [TTT e). we plot the evolution of the neutrino luminosity for models E 
and E2. In contrast to the results for a < 0.75M, the luminosities for C = 1 &iid 
1/3 differ significantly. The reasons are that (i) for a = 0.9M, the rotation radii of 
the torus can be small enough {uj ~ 2.5M; see Table II) to constitute a compact 
torus and that (ii) the accretion time scale is long enough to halt quick accretion for 
large-spin black holes. Consequently, the density can increase to ^ 10^^ g/cm'^ even 
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with M « O.25M0 [see Fig. [ToTb)]. Because the optically thick region is large for 
such accretion flow, the decrease of (" significantly increases the size of the optically 
thin region and the luminosity. 

In Fig. fTlT f). the evolution of the conversion efficiency, L^/M^,, is plotted for 
models E and E2. It is seen that the conversion efficiency for model E2 is about 
2.5 times as large as that for E. For a = 0.9M, the hypothetical maximum conver- 
sion efficiency is about 15.5% (cf. Table II). For model E2, the average conversion 
efficiency {AEu/AM:f,) is ~ 6%, which is thus ~ 40% of the maximum value. For 
a ~ M, approximately 40% of the rest-mass energy may be converted into thermal 
energy.^^) Our results for model E2 suggest that a conversion efficiency of ~ 20% 
may be possible for o ~ M. 

We determine the neutrino optical depth in a simple qualitative way, using 
the local density. In an actual system, the optical depth depends on the density 
distribution, temperature profile, neutrino energy, geometry of the black hole, and 
path of neutrinos in the curved geometry. Although the assumption = Cpii is not 
bad qualitatively, the error is not small quantitatively: Comparing the results for 
C = 1 and 1/3, we infer that the error in the luminosity is approximately a factor of 
2-3 for large values of the black hole spin. To obtain the neutrino luminosity more 
accurately, it is necessary to adopt a more sophisticated method for determining the 
optical depth, in particular for large values of a. This is beyond scope of this paper 
and left for future study. 

§5. Summary and discussion 

5.1. Summary 

We have reported our first numerical results of a GRMHD simulation for neutrino- 
cooled accretion tori. We solved the GRMHD equations in the fixed gravitational 
field of Kerr black holes of mass 4Mq with a realistic EOS and with neutrino cool- 
ing. The simulation was carried out systematically for a wide range of values of the 
black hole spin and the torus mass. Below we summarize the results of the numerical 
simulation. 

• In the presence of both poloidal magnetic fields and differential rotation, the 
magnetic field strength is amplified by the winding of the field lines and by the 
MRI until the electromagnetic energy reaches ~ 10% of the rotational kinetic 
energy. Then, the magnetic stress induces angular momentum transport via 
magnetic braking and the MRI, resulting in a quasistationary accretion onto 
the black hole. It also drives turbulent motion of matter, which subsequently 
generates shocks that convert kinetic energy into thermal energy. Through these 
processes, the temperature of the torus increases typically to ~ 10^^ K. The 
ratio of the electromagnetic energy to the rotational kinetic energy is maintained 
at ~ 10% in the quasistationary accretion phase. This electromagnetic energy 
is comparable to the internal energy in the quasistationary phase. 

• The maximum density of the torus takes values in a wide range, between ~ 
lO^'^ g/cm^ and ~ 10^^ g/cm^, depending on black hole spin a and mass of the 
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torus M^. For larger values of a, the maximum density tends to be higher for a 
given value of M^,, because (i) the location of the ISCO is closer to the horizon, 
leading to a more compact torus, and (ii) the accretion rate is suppressed by 
the small ISCO radius, halting the infall of the matter and resulting in the 
formation of a higher-density region near the ISCO. For a larger torus mass, 
the maximum density becomes larger. 

• In the case that the density is sufficiently high, some of the neutrinos are trapped 
in the accretion flow, which suppresses the neutrino emission rate. This tends 
to happen for large values of a and for large values of the torus mass. 

• Before the torus relaxes to a quasistationary state, the accretion rate reaches 
~ IOM0 /s, but after it relaxes, a typical accretion rate is Mq /s. The 
corresponding accretion time scale is ~ 50 ms in the early phase, but it relaxes 
to 100-200 ms in the quasistationary state. This accretion rate is in agreement 
with that in the case that an a- viscosity in the range = 0.01-0.1 is included. 

• The maximum neutrino luminosity is a few xlO^'^ crgs/s for the torus mass 

O.25M0, irrespective of the value of the black hole spin for < a/M < 0.9. 
(For a = 0.9M and ( = 1/3, the maximum reaches exceptionally 10^^ ergs/s.) 
In the quasistationary phase, it is between 10^^ ergs/s and 10^^ ergs/s, which 
depends strongly on the black hole spin. The efficiency of the conversion to 
neutrinos, L^/M^,, is between 1 and 10 %. This value depends on the black hole 
spin and the effect of neutrino trapping. The total emitted neutrino energy is 
2 X 10^^-2 X 10'^^ ergs for O.25M0. This also depends strongly on the 

black hole spin and the effect of neutrino trapping. 

• The neutrino luminosity, L^, and the total emitted neutrino energy depend 
strongly on the mass of the torus, M*. The maximum value of the luminosity 
and the total emitted energy are approximately proportional to M* for M* ^ 
OAMq. 

• Neutrinos are emitted efficiently in the region with w ^ lOM and \z\ ^ 5M. 
In particular, the neutrino emissivity is highest near the inner surface of the 
accreting torus. 

• The neutrino luminosity and the conversion efficiency for a ^ 0.75M are larger 
by a factor of ^ 2 than those for a = 0. However, moderate values of the spin 
(a ^ 0.5M) do not help to significantly increase the luminosity for a torus of 
mass = O.1-O.4M0. 

• If the accretion flow is optically thin with respect to neutrino transport, the 
efficiency of the conversion to neutrinos may be larger than ~ 10% for a ^ 0.9M. 

5.2. Implications for GRBs 

Because neutrinos are emitted most efficiently near the inner surface of accreting 
torus, neutrino-antineutrino pair annihilation is expected to occur in the vicinity of 
the rotation axis of black hole. As a result, a pair plasma of electrons and positrons 
may be generated, forming a fireball near the rotation axis. If the energy density of 
this fireball is sufficiently high, it can drive GRBs.^^ Setiawan et al.^-* estimated the 
luminosity of gamma rays, and found that it could be ~ 10^°-10^^ ergs/s for 
Lu ~ 10^^ ergs/s. Such a large value would be sufficiently high to generate GRBs if 
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the total baryon mass in the vicinity of the rotation axis is sufficiently small. We 
find that the neutrino luminosity is ~ 10^^ ergs/s for a torus mass of ^ 0.2Mq and 
black hole mass of 4M0, irrespective of the black hole spin. It is believed that such a 
system may act as the central engine of GRBs (specifically, GRBs of short duration). 
In particular, for large values of a close to unity, the efficiency of the conversion to 
neutrinos could be high, and a strong GRB would then be driven. 

In our present treatment, the heating of matter by neutrinos emitted from the 
torus and ncutrino-antincutrino pair annihilation are not taken into account. In an 
actual system, neutrino heating will help the formation of matter outflow, which will 
contribute to sweeping baryons around the rotation axis. 

The neutrino luminosity decreases as a function of the torus mass. For a torus 
of mass ~ IO^^Mq, the luminosity is likely to be at most ~ 10^^ crgs/s. Because 
the pair annihilation rate of a neutrino-antineutrino process is proportional to the 
square of the neutrino luminosity, Ei,p is likely to be ^ 10^^ ergs/s, and hence, tori 
of small mass arc imlikcly to drive the observed GRBs. Numerical simulations have 
shown that the merger of binary neutron stars of sufficiently large mass can produce 
a system consisting of a rotating black hole and a torus of spin a/M = 0.7-0.8 
(e.g., Refs. 45), 46)). However, the torus mass is in general small, <^ O.IMq, unless 
the mass ratio of the two neutron stars is sufficiently small, ^ 0.7-0.8. Even for a 
small mass ratio of ~ 2/3, the torus mass is at most ~ O.IMq. This indicates that 
the merger of binary neutron stars with sufficiently large total mass could result 
only in weak GRBs in most cases. If the total mass of binary neutron stars is not 
large enough for the direct formation of a black hole, a hypermassive neutron star is 
formed. ^^-^ (See Ref. 47) for the definition of the hypermassive neutron star.) If it is 
strongly magnetized, such a hypermassive neutron star will collapse and become a 
rotating black hole with a massive torus, due to the transport of angular momentum 
induced by magnetic effects. ""^^^ It is believed that the resulting system of a black 
hole and a torus may be the central engine of GRBs, in contrast to the case that 
the total mass of the binary neutron stars is large enough for the direct formation of 
black hole. Another possibility as the source of GRBs is the merger of black holes 
and neutron stars. If the mass of the black hole is not large (^ 4Mq), the resulting 
mass of the torus around the black hole can be ^ O.IMq. ^^^^^^^ Therefore, it is 
regarded as a strong possibility as the source of GRBs. 

A characteristic feature found in the MHD simulation is that the luminosity 
varies on a short time scale of a few ms, in contrast to the results reported in Refs. 
8) and 9). This is due to the fact that shock heating in the magnetized accretion 
torus is associated with turbulent matter motion driven by a highly variable magnetic 
stress. The light curve of GRBs is often highly variable. If the GRBs are driven by 
pair annihilation of neutrino-antineutrino pairs emitted from the magnetized torus, 
the luminosity curve of neutrinos should also be highly variable. Thus, the light curve 
of GRBs may be naturally explained if the central engine is composed of a rotating 
black hole and a magnetized accretion torus of mass ^ O.IMq. 
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5.3. Future tasks 

This work is the first step toward obtaining a detailed understanding of the evo- 
lution of dense, hot magnetized tori surrounding a stellar-mass black hole. There 
are several tasks left for the future. One is to take into account more realistic mi- 
crophysics, incorporating more sophisticated EOSs. To this time, two tabulated 
finite-temperature EOSs have been published for public use. These incorporate the 
effects of heavy nuclei as well as free nucleons and a-particles.^''-*'^'^^ In these EOSs, 
the pressure, internal energy, and other various quantities are written as functions 
of the density, electron fraction, and temperature, like the EOS used in this paper. 
Thus, it would be straightforward to change the numerical code to use more realistic 
tabulated EOSs. Numerical results, such as the rest-mass accretion rate, temper- 
ature, and neutrino luminosity, may depend on the chosen EOSs. To clarify this 
dependence, we plan to perform simulations with such EOSs. 

Another task is to improve the scheme for treating neutrino transfer. As pointed 
out in Refs. 2)-4),8) and 9) and reconfirmed in this paper, the luminosity is sig- 
nificantly affected by neutrino-trapping effects. In this paper, we determined the 
neutrino optical depth by simply using the local density. However, the actual optical 
depth depends on the density distribution, temperature profile, neutrino energy, and 
geometry of the black hole. As shown in mA\ the luminosity depends strongly on the 
chosen optical depth for large values of a and M*, because such tori are composed 
both of optically thick and thin regions. To provide a better estimate of the neutrino 
luminosity, we have to treat the optical depth more carefully. 

Elucidating the exact neutrino trajectories in the curved spacetime is also an 
important problem. Neutrinos are most efficiently emitted in the vicinity of the 
black hole. This implies that a non-negligible fraction is swallowed by the black 
hole due to large curvature effects. We plan to study this effect in the future. An 
associated problem is to estimate the annihilation rate of neutrino and antineutrino 
pairs. As mentioned above, this process generates a pair plasma for GRBs. Although 
Setiawan et al. already calculated this rate in their simulation,^) they did not take 
into account the curvature effect. We consider that estimating the energy generation 
rate in curved spacetime is an important problem. 

The final problem that warrants further study concerns nonaxisymmetric ef- 
fects. The development of the MRI in an axisymmetric system could be different 
from that in the nonaxisymmetric, 3D case.^^^ Turbulence tends to persist more in 
the 3D case due to the lack of symmetry. Specifically, according to the axisymmetric 
anti-dynamo theorem,^'''' the sustained growth of magnetic field energy is not possi- 
ble through axisymmetric turbulence, as demonstrated by numerical simulations.^^) 
In our present simulation, we found that the neutrino luminosity decreases with a 
time scale of ~ 100 ms, implying that the efficiency of shock heating in the torus 
decreases with time. This may be partly due the anti-dynamo effect. McKinney 
and Gammie^^) have performed axisymmetric simulations of magnetized tori accret- 
ing onto Kerr black holes and have found good quantitative agreement with the 3D 
results of De Villiers and Hawley^^^ for the global quantities E/M^ and j/M^ for 
t < 2000M. Thus, the results presented here likely provide (at least) a good qualita- 
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tive picture at least for the short-term evolution. However, to clarify the longer-term 
evolution, 3D simulations will eventually be necessary. 
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Appendix A 

Derivation of the equation for T^^ 

The total energy-momentum tensor, T^, is defined by the sum of energy mo- 
menta of all the matter fields, and obeys the conservation equation 

= 0. (A-1) 

In the context of the present paper, is composed of the energy momenta of 
baryons, electrons, positrons, radiation, and neutrinos. Then, we split the total 
energy-momentum tensor as 

where T^i, is written by Eq. (12-91) which is composed of baryons, electrons, positrons, 
radiation, and thermal neutrinos. In this paper, the energy-momentum tensor of 
thermal neutrinos are simply written by 

pe^u^Uy + Pyg^y. (A-3) 

On the other hand, is the contribution from free-streaming neutrinos for which 
the evolution equation is assumed to be 

v^rN't = Qu. (A-4) 

Then, the equation for T^^y is 

V;.n = -Qu. (A-5) 

In the present simulation, in which the MHD equations are solve in the fixed 
background of a black hole, we do not have to determine T^, because such a term 
never appears in the basic equations, and thus, we do not solve Eq. (|A-4p in this 
paper. In fully general relativistic simulation in which the Einstein equation is solved, 
however, it is necessary to solve this equation, because it appears as the source terms 
in the equations of the geometric variables. 

Appendix B 

Derivation of equation for electron capture rate 

In this appendix, we derive the equations for the electron capture rate (|2-56p 
and the associated energy emissivity (|2-60p . 
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The electron capture rate by nuclei, Ag-, is in general given by [e.g., Ref. 53)] 
_ln2^ (2J. + l)e-^-/(^^) ^ 

where the sums over i and j run over states in the parent and daughter nuclei, 
respectively. The constant K (= 6146 ± 6[s]) is defined as 

K= ^'^'ff 1 . (B.2) 

where Gf is the Fermi coupling constant, V^d is the up-down element in the Cabibbo- 
Kobayashi-Masukawa quark-mixing matrix, and i^v = 1 is the weak vector coupling 
constant. Further, the quantity G{Z, A,T) = exp(— S,/A;T) is the partition 
function of the parent nucleus, Mij is the reduced transition probability, and fij is 
the phase space integral, given by 

^•=(^^j J x\x + Qj?F{Z,x)S,-{x)[l-S,{x + Cij)]dx. (B-3) 

Here x is total energy of electrons divided by kT, and Qj is the difference between 
the nuclear mass-energies of the ground states of the parent and daughter nuclei in 
units of kT: 

Cij = {Mp^rc^ - M^^y + Ei-Ej). (B-4) 

In this expression, Mpar and M^au are the nuclear masses of the parent and daughter 
nuclei, respectively, and Ei and Ej are the excitation energies of the initial and 
final states. The lower limit of the integral, rji, is the capture-threshold total energy 
in units of kT, which is given by rji = uieC^/kT if Qj + mf.(? > and rji = \C,ij\ 
otherwise. Also, in ()B-3p . Se and Si, are the electron and positron distribution 
functions, which are assumed to be described by the Fermi-Dirac distributions with 
(matter) temperature T and chemical potential i]^ = j-i/kT, 

Se±{x) = -—i , Su{x) = ^ i . , (B-5) 

"*?e± -)- 1 -|- 1 



X- 



where ry^ and rj^ are the electron (positron) and neutrino chemical potentials in 
units of kT. The quantity F{Z, x) is the Fermi function, which corrects the phase 
space integral for the Coulomb distortion of the electron wave function near the 
nuclei. 

For actual calculations of the capture rate, we follow a prescription introduced by 
Fuller, Fowler, and Newman. ^^^'^^^ In this prescription, the capture rate is rewritten 
in terms of the effective f-t value {ft) as 

Ae-=ln2^, (B-6) 
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where we have factored out the contribution of the Fermi function from the phase 
space integral [i.e., lij = fij/{F(Z,x))] and have put this contribution into the 
effective f-t value. The reduced phase space factor is then 



kT 
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+ 1 



dx. 



(B-7) 



In terms of the relativistic Fermi integrals Fk{r]), the phase space factor can be 
written 

kT \^ 1 



m^c^ 



1 



-le, 



(B-8) 



where 



/•oo 

Ie= x'^{x + Cij) 



+ (20, +4,?^) [Fsiv!- 



dx 



F2{€ - Q: 



+ [(r?^)^ + 2ir^^fQ, + [v'^fiCfjf] [Fo{C - - Fo{r^^ - - t?^)] . (6-9) 

For electron capture by free nucleons {p + e~ n-\- v^, ri + e"*" p + i>e), the 

effective f-t value is {ft) = 10^°^^ s,^^) and the difference between nuclear mass- 
energies is 

Qj = - rrinc^ + Ep- En). (B-10) 

In this paper we have assumed that the nuclear mass-energy difference as well as elec- 
tron mass is much smaller than the total energy of relativistic electrons. Accordingly, 
the phase space factor is simplified to 
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In the region where neutrinos freely stream out, the phase space factor is further 
simplified to 

^V^4(r,f-). 



•■tj 



meC'' 



(B-12) 



The electron capture rate by free protons is given by 



K- = Kc 



Similarly, the positron capture rate is 
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(B-14) 
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The energy emission rates of neutrinos in units of nieC^ associated with the 
electron capture by free nucleons are given by 



where 



kT 
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Similar calculations to derive equations ()B-13P and ()B-14p give 
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dx. (B-16) 



(B-17) 



To avoid the direct numerical integration of the relativistic Fermi integrals Fj (rj) , 
we adopt approximate expressions for them following Refs. 56) and 55). For j = 0, 
the relativistic Fermi integrals can be integrated exactly, and we obtain i^o(^) = 
ln(l + e'') = r] + ln(l + e~^), which immediately gives 



Then, using the well-known recursion relations 

Fj{7]) = Fj{0)+j r Fj-i{x)dx, 
Jo 

F_,(r?) = Fj{0)-j r Fj^ii-x)dx, 
Jo 



(B-18) 



(B.19) 



and noting that Fj{0) = (j!)(l — 2 ■^)C{j + l), where ^{x) is the zeta function, we 
find 
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The standard expansion of the relativistic Fermi integrals for < is given by (e.g., 
Ref. 57)) 



1=0 



(B-21) 



We adopt the first three terms for our approximation of the relativistic Fermi inte- 
grals for ?7 < 0: 



F^irj) « 24 
F^iv) « 120 
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(B-22) 
(B-23) 
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Our approximation of the relativistic Fermi integrals for 77 > is obtained as follows 
by using relations (|B-20p and imposing continuity conditions on the values and their 
first derivatives: 

F4(r?) w 45.597/ + — t?^ + -t + ^4(-??) 

6 5 

F5 iv) ^ 236.65 + —7j^ + —r]^ + - F^i-T]). (B-24) 

ODD 
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